/// Script to produce relative risk graph and estimates using nlcom command and coefplot package
/// Fabien Cottier
/// SPEI DATA


///////// Graph immediate SPEI effect only(SPEI Y0)

* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est	
	local models ""

	* Loop to create multiple nlcom statement
	forvalues j=-.8(0.1).8 {

		local x=round(`j',.1)

		if "`2'"=="linear"{
			local models="`models' (_b[wmean_speiq]*`x'+_b[L1.wmean_speiq]*`x'+_b[L2.wmean_speiq]*`x'+_b[L3.wmean_speiq]*`x')"
		}
		
		else if "`2'"=="quadratic"{
			local models="`models' (_b[wmean_speiq]*`x'+_b[wmean_speiq#wmean_speiq]*(`x')^2+_b[L1.wmean_speiq]*`x'+_b[L1.wmean_speiq#L1.wmean_speiq]*(`x')^2+" ///
				+"_b[L2.wmean_speiq]*`x'+_b[L2.wmean_speiq#L2.wmean_speiq]*(`x')^2+_b[L3.wmean_speiq]*`x'+_b[L3.wmean_speiq#L3.wmean_speiq]*(`x')^2)"
		}
		
	}

	* estimates relative risk using nlcom
	est resto `1' 
	local df_est=e(df_r)
	nlcom `models', post df(`df_est') 
	estimates store est_rr
	
	* produce plot using coefplot
	coefplot (est_rr), ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+ 75%") xlabel(1.5 "-0.75" 4 "-.5" 6.5 "-0.25" 9 "0" 11.5 "0.5" 14 ".5" 16.5 ".75" ) ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ciopts(recast(rarea)) yline(1) 
end




///////// Split sample graph effects SPEI (SPEI Y0)


* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_Sp	// MuL_Sp for multiple lag with split sample
	local models_L0 ""  

	* Loop to create multiple nlcom statement
	forvalues j=-.8(0.1).8 {

		local x=round(`j',.1)

		if "`3'"=="linear"{
			local models_L0="`models_L0' (_b[wmean_speiq]*`x'+_b[L1.wmean_speiq]*`x'+_b[L2.wmean_speiq]*`x'+_b[L3.wmean_speiq]*`x')"
		}
	
		else if "`3'"=="quadratic"{
			local models_L0="`models_L0' (_b[wmean_speiq]*`x'+_b[wmean_speiq#wmean_speiq]*(`x')^2+_b[L1.wmean_speiq]*`x'+_b[L1.wmean_speiq#L1.wmean_speiq]*(`x')^2+" ///
				+"_b[L2.wmean_speiq]*`x'+_b[L2.wmean_speiq#L2.wmean_speiq]*(`x')^2+_b[L3.wmean_speiq]*`x'+_b[L3.wmean_speiq#L3.wmean_speiq]*(`x')^2)"
		}
		
	}
	
	* estimates relative risk using nlcom // With High labor agriculture
	est resto `1'
	local df_est=e(df_r)
	nlcom `models_L0', post df(`df_est')
	estimates store est_rr_L0H

		
	* estimates relative risk using nlcom // With Low labor agriculture
	est resto `2'
	local df_est2=e(df_r)
	nlcom `models_L0', post df(`df_est2')
	estimates store est_rr_L0L
	
	* graph
	coefplot (est_rr_L0H), bylabel(High Agriculture - Immediate effect) || ///
		(est_rr_L0L), bylabel(Low Agriculture - Immediate effect) ||, ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+ 75%") xlabel(1.5 "-0.75" 4 "-.5" 6.5 "-0.25" 9 "0" 11.5 "0.5" 14 ".5" 16.5 ".75" ) ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ///
		ciopts(recast(rarea)) yline(1)  byopts(row(3))
		
end


